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Abstract 

We numerically compare the semiclassical "frozen Gaussian" Herman-Kluk propagator [Chem. 
Phys. 91, 27 (1984)] and the "thawed Gaussian" propagator put forward recently by Baranger 
et al. [J. Phys. A 34, 7227 (2001)] by studying the quantum dynamics in some nonlinear one- 
dimensional potentials. The reasons for the lack of long time accuracy and norm conservation 
in the latter method are uncovered. We amend the thawed Gaussian propagator with a global 
harmonic approximation for the stability of the trajectories and demonstrate that this revised 
propagator is a true alternative to the Herman-Kluk propagator with similar accuracy. 
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I. INTRODUCTION 



The inclusion of quantum effects into classical dynamical calculations is of fundamental 
as well as of practical interest. It may, e. g., help in understanding the nature of interference 
related effects and, for systems with many degrees of freedom (DOF), it may make quantum 
calculations possible in the first place. The most fundamental approach to reach this goal 
is based on the application of the semiclassical approximation to the quantum mechanical 
path integral expression for the time-dependent propagator A major drawback of the 
standard time-dependent semiclassical Van Vleck-Gutzwiller formula P, 0] is the so-called 
root search problem. One has to find solutions to classical boundary value problems, i. e. 
classical trajectories connecting a given initial with a given final position in a specified time. 
This may become numerically cumbersome for longer times, especially if many DOF are 
involved in the classical dynamics, as is usually the case. 

Initial value representations (IVRs) of the semiclassical propagator have therefore been 
the focus of interest for quite some time . In the IVRs, root searches are 

avoided by transforming the semiclassical propagator into an integral expression leading to 
the emergence of classical initial value problems instead of boundary ones. Specifically the 
He.n,an-Kluk (HK) IVR e.p.e=.on based on fl.ed wHth, . e. "frozen" G— has been 
quite successfully applied since its first derivation Recently, an alternative formulation, 
based on "thawed" Gaussians, has been brought back into the limelight in Refs. 11 1 an d |l2l| . 
The single trajectory version, for a propagated Gaussian in the formulation of Ref. [uj is 
equivalent to Heller's Gaussian wave packet dynamics 3]; evolving a Gaussian wave packet 
by using only its center trajectory and allowing the width of the Gaussian to "breathe". The 
full integral expression based on thawed Gaussians (i. e. the multitrajectory version, see also 



Ref. 



12| |) has already been discussed by Kay ,8j and has been labelled thawed Gaussian 



approximation (TGA). Rather surprisingly, Kay found first numerical evidence that the 
TGA is numerically far inferior as compared to the HK frozen Gaussian Approximation (HK 
FGA). As the main result of this article, we propose a "breathing", individually determined 
from the global properties of the classical dynamics on the energy shell which is defined for 
each trajectory by its initial point in phase space, in order to adjust the shortcomings of 
TGA. 

In the following, we will explicitly examine several examples to shed some light on the 
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reason for the failure of the original TGA. Simultaneously, it will become clearer why FGA 
is numerically superior. To set the stage, in Sec. II, we briefly review Kay's results [? ] of a 
general semiclassical IVR expression based on coherent states. In Sec. Ill, the performance 
of the two special cases, HK FGA and TGA, is contrasted for two different anharmonic test 
potentials. In Sec. IV, the reason for the TGA failure is underlined by a simple improvement 
before we introduce the global harmonic approximation which improves TGA to the level 
of FGA. Conclusions and an outlook are given in Sec. V. 



II. SEMICLASSICAL INITIAL VALUE REPRESENTATIONS 

A general expression for the quantum propagator, subsuming many of the previously 
known semiclassical initial value representations, has been derived via the principle of sta- 
tionary phase equivalence by Kay 8]. He used the fact that the Ansatz as an integral over 
phase space has to lead to the Van Vleck-Gutzwiller form, if the integration is performed in 
the stationary phase approximation. 

Several early semiclassical approaches, e. g. the ones of Heller and Herman and Kluk 
are based on Gaussian wave packets of the form 



in coordinate representation. In his seminal work 8], Kay has used these states with arbi- 
trarily time-dependent width parameters 7(t) in order to construct a general semiclassical 
integral expression for the propagator. For our purposes, we need only the special version 
without a cross term between initial and final Gaussian, given by 

K{x,t;x',0) = / -^^-^{x\g^,{pt,qt))Rip^,qi,t)exp{iS{pi,q^,t)/h}{g^^{pi,qi)\x'), (2) 



where the pre-exponential factor 



1 f dqt ^ dpt dqt 1 dpt 

717^ + 72 7^ ^"-7172 7^ 



1/2 



(3) 



_2A/Re7iRe72 \ % dpi dpi ih dq, 

is determined by the stability matrix (or monodromy matrix) elements defined in the matrix 
equation 

^qt \ _l dqt/dqi dqt/dpi \ / 6q, 
6pt / I dpt/dqi dpt/dpi I \ 6p., 



(4) 



which connect the initial with the final deviations of trajectories {pt, qt) in phase space, while 
the phase is determined by the classical action S{pi,qi,t) = Jq Ldt' being the time integral 
over the Lagrange function. 

From the general expression, Eq. the HK approximation is recovered by setting 
7i = 72 = 7 with a real positive constant 7. The explicit form of the HK FGA prefactor is 
then 



"^2 \dqi dpi dpi ih'y % / ' 



while by setting 



_ i dpt/dqi + ihjdpt/dpi , , 

r—,^'^ h dqt/dqi + ihjdqt/dpi 
again with real positive 7 jlj] and 72 = 7, we recover the TGA. This semiclassical IVR has 

recently also been discussed by Baranger et al. jll|, [3| • The prefactor turns out to be 

In the following, we will numerically investigate the TGA expression and compare the 
results to FGA and full quantum results. Later on, the general expression ^ will then also 
allow us to search for improvements of the TGA propagator. 



III. LONG TIME ACCURACY OF HK FGA VERSUS TGA 

For reasons of simplicity, the wave function to be propagated is chosen to be a Gaussian 
\l/(x,0) = {x\g^{po,qo)) localized around {pQ,qo) in phase space and having the same width 
)arameter as the initial Gaussian in the integral expression of the propagator in Eq. 
3]. In the following we will compare the overlap between the time propagated and the 
unpropagated wavepacket, the so-called autocorrelation function 

/oo 
c/x^*(x,0)^(x,t), (8) 
-00 

for the different semiclassical approaches. 

Inserting the propagator from Eq. (0), we obtain the semiclassical autocorrelation func- 
tion 

c{t) = f ^^^^{gj{po,qo)\9^i{Pt,qt))R{Phqi,t) 

exp{iS{pi,qi,t)/h} {9^2{Ph<li)\9^{Po,qo))- (9) 



The overlap integral between the Gaussians can be performed analytically, whereas the phase 
space integration in Eq. © is done numerically by the Monte Carlo integration procedure 
described in Ref. To test the numerical performance of the different methods we chose 
two different nonlinear potentials which we now introduce. 

First, we employ the Hamiltonian taken by Baranger et al. |1JL,] to test their integral 
kernel, 
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Hb = ^ + 2^06""^ cosh(aa;). (10) 

Throughout the rest of the article, we use only this classical version of the Hamiltonian 
(without Gaussian smoothing) with (scaled) units, particularly scaled time ts, indirectly 
defined by the parameters taken from Ref. jn| as/x = l,Vo = l,A = 5,Q; = l, and a 
Planck constant oih = 0.05. The initial conditions of the wave packet are go = and = 1 
and 7 = 100/9. "Smoothing" of the Hamiltonian by using its coherent state matrix element 
leads only to minor changes (not shown) in the results to be presented below. 

As the second test bed, we choose the well known Morse potential with the Hamiltonian 

HM = Y^ + Vo{l-e-^^f (11) 

which is a reasonable model for a chemical bond. The parameters taken from Ref. j3| as 
A = 0.08, Vo = 30, which together with yU = 1 define a scaled time ^m- As initial conditions 
for the wavepacket we define Q'o = Po = and 7 = 12. 

In Fig. we compare different semiclassical IVR results (i. e. using the full integral 
expression) to the quantum result, calculated by the split operator method, for the two model 
potentials. For the Baranger potential, we have calculated the autocorrelation function up 
to scaled times of 110, as displayed in Fig.^a). In accord with the decay of the norm in Fig. 
^c), also the autocorrelation function shows a strong decay, not present in the full quantum 
curve. We state explicitely that the decay of the norm for TGA and also the discrepancy 
of the norm from unity in the HK FGA case are not in contradiction to the fact that both 

n 

methods are unitary in the (analytical) stationary phase sense (see e. g., Ref. |17[ for the 
HK FGA case). 

Nevertheless, for the Baranger potential, the quality of the numerical TGA results is 
disappointing. For the very short times (up to 10 in scaled units) investigated in Ref. 
h|, even the one trajectory (kernel) version of TGA, i. e.. Heller's Gaussian wave packet 
dynamics, does better than the full integral version (not shown). This method is, however. 



not suited for longer times due to the nonlinearity of the potential. In contrast, as can also 
be observed in Figs, ^a) and (c), the HK FGA approach gives a reasonable agreement 
with the full quantum solution. 

The same kind of discrepancies of TGA from the quantum split operator calculation 
as for the Baranger potential can be observed in the case of the Morse potential. The 
corresponding results are given in Figs. [T] (b) and[T] (d). As has been noted already in jl^ . 
for the Morse potential, the HK FGA approximation represents the fine details of the full 
quantum result very well. The multitrajectory TGA expression, however, already fails again 
at intermediate times. 

IV. IMPROVEMENTS OF THE TGA 

With the comparison of the different results at hand, we have the means to find the 
reason for the poor performance of the TGA expression. 

A. Generalized TGA 

The reason for the disappointing TGA results becomes clear by taking a closer look at 
the time-dependent width parameter 71. To this end, we concentrate on the Hamiltonian 
of Eq. ()10|) with the Baranger potential. As can be seen in the lower panel of Fig. |21 the 
real part of the complex width parameter is almost always close to zero for longer times. 
This leads to a strongly oscillating exponent in the TGA integral expression. Although the 
Monte Carlo results which we present are converged, the strong oscillation of the integrand 
leads to the rapid deterioration of the results. The reason for the satisfactory behavior of 
the HK FGA in the cases we investigated is obviously the fixed real width parameter which 
by definition does not decay. 

This observation opens up a possibility to improve on the standard TGA expression in 
a simple way. If one could render the width parameter's "destructive influence" on the 
results less dramatic, this would certainly ameliorate the results. An easy way to achieve 
this goal is by taking (higher-order) roots of the width parameter 71. This is done in the 
upper two panels of Fig. |21 where in panel (b)-the square root -and in panel (a) the fourth 
root of the width parameter in panel (c) is taken. Clearly, the decay of the real part is 
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less dramatic. The effect that this change has on the autocorrelation function can now be 
checked by inserting the new width parameter into the general integral expression of Eq. ^ 
and calculating the resulting prefactor. The emerging propagator is then still semiclassical 
in the sense of stationary phase equivalence. The numerical results are shown in the upper 
panels of Fig. |21 Already in the case of the fourth root (see panel Fig. 121(a)), the generahzed 
TGA result is almost as good as the HK FGA. By making the final Gaussians narrower in 
width, we have reached an intermediate goal. TGA can now better capture the real quantum 
dynamics. 

There is, however, one further aspect to this improvement. The issue whether the gen- 
eralized TGA can do as well as HK FGA with fewer trajectories. To study this, we have 
investigated the convergence of the Monte Carlo integration for FGA, the two different gen- 
eralized TGAs, and the standard TGA. For the Baranger potential example, all methods 
require on the order of several times 10^ trajectories to achieve convergence. The HK FGA 
needs a slightly higher number of trajectories than standard TGA. Taking roots of higher 
and higher order makes the width almost constant and, therefore, the TGA propagator ap- 
proaches the FGA. This also leads to more trajectories which are required for convergence. 

Besides taking roots of the width parameter, there is yet another way to improve the TGA, 
making it a true alternative to the FGA. This approach will be studied in the remainder of 
this article. 



B. Global harmonic approximation 

A starting point for an improvement of the TGA, different in spirit from the one presented 
above, is based on the fact that all time- dependent semiclassical approaches become exact 
for the case of the harmonic oscillator. The TGA is especially favorable for the case of a 
harmonic potential, whose frequency u does not coincide with the width parameter of the 
Gaussians used in the integral expression for the propagator. Already the single trajectory 
Gaussian wave packet dynamics gives the exact analytical solution (the so-called squeezed 
state) which breathes (for uj ^ ^"^1 1^) according to the time-dependent width parameter 

i — fiu sm{ut) + ih'j cos{ujt) 

~ • (-'-^J 

h cos{ujt) + ih^y/i^fiuj) sm{ut) 
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The individual Gaussians of TGA breathe in the same way as the solution does, allowing 
the convergence of the numerical results to be achieved with fewer trajectories than in the 
HK FGA case, also explaining the remarks at the end of the previous subsection. Is this 
observation of any use for realistic non-harmonic potentials? 

The answer is yes and we want to demonstrate this for the case of the Morse oscillator 
with the parameters given in the previous section. For the bound motion in any potential 
we can define an action variable 



by integrating over a closed contour C in phase space. The Hamiltonian does not depend 
on the corresponding angle variable and can be expressed as a function of /. Its derivative, 
with respect to the action variable, then defines a frequency 

which in anharmonic potentials depends on energy and, in the simple one DOF case, can be 
determined analytically. The idea of the improvement of TGA is now to use the oscillation 
frequency given above as the one of a "fictitious" harmonic oscillator and let every Gaussian 
breathe with the corresponding 71. While all stability matrix elements for the calculation of 
the semiclassical propagator are taken to be the harmonic ones, the trajectories themselves 
are propagated using the full nonlinear force field. For every trajectory, the approximation 
is a global harmonic approximation, with respect to the stability information. Note that the 
present approach is different from the the so-called nonlinear Gaussian wave packet dynam- 
ics, introduced by Tomsovic and Heller which uses a local expansion of the potential 
up to second order (as Heller's original Gaussian wave packet dynamics) and is not an IVR 
formulation but a root search method in spirit. 

We are now using the global harmonic approximation for the width parameter and also 
the prefactor in Eq. Q to calculate the autocorrelation function in the Morse oscillator case. 
The results are shown in Fig. EJ A dramatic improvement toward the HK FGA result can be 
observed. The reason is that now 71 behaves as in the harmonic case, i. e. it does not decay 
to zero for later times, thus leading to reasonable results. We found that the same kind 
of improvement toward HK is also present for wave packets with different initial conditions 
(not shown). 




(13) 
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V. CONCLUSIONS 



Based on numerical studies for two different anfiarmonic potentials (Baranger and Morse), 
we have identified the decay of the width parameters of the Gaussians in time as the reason 
for the failure of the TGA to provide a good semiclassical propagator. Following this notion, 
we have introduced the global harmonic approximation which mimics for each trajectory 
its stability properties in time according to those of a harmonic oscillator. Its frequency is 
defined through an action-angle analysis of the dynamics on the energy shell set by the initial 
conditions of each trajectory. The TGA, amended by the global harmonic approximation 
for the stability of trajectories, is comparable in accuracy to the HK FGA propagator. 

More detailed studies, e. g. of the convergence properties of the different methods, on the 
choice of an optimal width parameter 7 for FGA, and on the applicability of the improved 
TGA also to unbound problems are worthwhile topics for future investigations. 
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FIG. 1: Comparison of different semiclassical versus quantum autocorrelation functions and 
tlie semiclassical results for the norm for the Baranger potential (a) and (c) and the Morse 
potential (b) and (d) with the potential parameters given in the text. Full line: quantum 
result, dashed line: HK FGA, dotted line: TGA 

FIG. 2: The real part of the time dependent width parameter in the TGA expression: (c) 
standard 71, (b) 71 = ^7!, (a) 71 = 71''^. 

FIG. 3: The standard and generalized TGA results in the case of the Baranger potential 
according to the prefactors in Fig. |21 Full line: quantum result, dotted line: TGA 

FIG. 4: The global harmonic TGA result compared to the HK FGA results for the Morse 
oscillator. Dashed line: HK FGA, dotted line TGA, and full line: global harmonic TGA 
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FIG. 1: Comparison of different semiclassical versus quantum autocorrelation functions and the 
semiclassical results for the norm for the Baranger potential (a) and (c) and the Morse potential 
(b) and (d) with the potential parameters given in the text. Full line: quantum result, dashed line: 
HK FGA, dotted line: TGA 
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FIG. 3: The standard and generalized TGA results in the case of the Baranger potential according 
to the prefactors in Fig. |2 Full line: quantum result, dotted line: TGA 
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